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1 Introduction 

Molecular dynamics (MD) is the ultimate method to gain detailed atomistic information of dy- 
namical processes that are difficult to access experimentally. However, an important bottleneck 
of atomistic simulations is the limited system- and timescales. Depending on the complexity 
of the forcefields (Ab Initio MD being extremely more expensive than classical MD) systems 
typically consist of 100 to 100000 molecules that can be simulated for a period of nanoseconds 
till microseconds. Therefore, many activated processes can not be studied using brute-force 
MD because the probability to observe a reactive event within reasonable CPU time is basi- 
cally zero. Typical examples are protein folding, conformational changes of molecules, cluster 
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isomerizations, chemical reactions, diffusion in solids, ion permeation through membranes, 
enzymatic reactions, docking, nucleation, DNA denaturation, and other types of phase tran- 
sitions. If these processes are treated with straightforward MD, the simulations will endlessly 
remain in the reactant states. Still, if an event would happen, it can go very fast. The time it 
actually takes to cross the barrier is usually much shorter than this computational accessible 
timescale. Therefore, rare event algorithms aim to avoid the superfluous exploration of the 
reactant state and to enhance the occurrence of reactive events. The methods that I will 
discuss are the reactive flux (RF) method [Tl[2] and the more recent algorithms that originate 
from the transition path sampling (TPS) [3HZ] methodology. These comprise the transition 
interface sampling (TIS) |8j and the replica exchange TIS (RETIS) [9l[T0], which are succes- 
sive improvements on the way reaction rates were determined in the original TPS algorithm. 
Partial path TIS (PPTIS) is an approximative approach in order to reduce the simulated 
path length for the case of diffusive barrier crossings. PPTIS is similar to Milestoning [12j . 
that was developed simultaneously and independently from PPTIS. For non-equilibrium sys- 
tems, the Forward Flux Sampling (FFS) was designed [13j. This method is based on the 
TIS formalism, but does not require prior knowledge on the phasepoint density. All these 
methods have in common that they aim to simulate true molecular dynamics trajectories at 
a much faster rate than naive brute force molecular dynamics. I will discuss the advantages 
and disadvantages of the different methodologies and introduce a few new relations and de- 
rive some known relations using a nonstandard approach. The descriptions of these methods 
given here are far from complete and, therefore, to obtain a more complete picture of the path 
sampling techniques I would like to recommend some very recent complementary reviews on 
these methodologies, [14H17j . In the end, I compare all the methods by applying them on a 
simple, though tricky, test system. The outcome illustrates some important pitfalls for the 
non-equilibrium methods that have no easy solution and show that caution is necessary when 
interpreting their results. 

2 Reactive Flux Method 

Low dimensional systems, such as chemical reactions in the gas phase, are usually well de- 
scribed by Transition State Theory (TST). TST assumes that the transitions from reactant 
to product state always follow a path on the potential energy surface such that it passes the 
barrier nearby the transition state (TS). The TS refers to the point on the energy barrier 
having to the lowest possible potential energy difference, with respect to the reactant state, 
that any trajectory must overcome in order to reach the product state. In this description, TS 
corresponds to a unstable stationary point on the potential energy surface having one imag- 
inary frequency (saddle point). In condensed systems, the saddle-point of the full potential 
energy surface is usually less meaningful. For instance, if we would consider the dissociation 
of NaCl in water, the TS would correspond to a state where the inter-ion distance is fixed to 
a critical value while all surrounding water molecules are full frozen into an icy state. It is 
needless to say that this does not correspond to our daily experience when dissolving some 
pinch of salt in a glass of water. TST can be generalized for higher dimensional systems using 
the free energy instead of the potential energy. The TS is then no longer a single point, but 
a multidimensional surface. In this case, the TST equation is determined by the free energy 
difference between the TS dividing surface and the reactant state. An important limitation of 
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TST is that the free energy barrier depends on the selected degrees of freedom that are used to 
describe the free energy surface. In addition, TST imphcitly assumes that any trajectory wih 
cross the free energy barrier only once when going from reactant to product state. Kramers' 
theory [TB] provides an elegant and insightful approach to correct for correlated recrossings if 
these originate from the diffusive character of the dynamics. However, there are several other 
sources for recrossings. For instance, if the selected degrees of freedom are not well chosen, the 
barriers in the free energy landscape do not always correspond the barriers in the underlying 
potential energy surface which ultimately determine the dynamics. The reactive flux (RF) 
method is able to correct for recrossings regardless their origin and is very powerful when 
when TST is close but not sufficiently accurate. 

The theory of the method originated from the early 1930s, far before the first applications 
of computers for molecular dynamics simulations [19]. After Wigner and Eyring introduce 
the concept of the TS and the TST approximation |20[l21j. Keck [22] demonstrated how to 
calculate the dynamical correction, the transmission coefficient. This work has later been 
extended by Bennett |23j. Chandler [24j and others [251[2B], resulting in a two-step approach. 
First the free energy as function of a single reaction coordinate (RC) is determined. This can 
be done by e.g. umbrella sampling (US) [27] or thermodynamic integration (TI) [2H]. Then, 
the maximum of this free energy profile defines the approximate TS dividing surface and the 
transmission coefficient can be calculated by releasing dynamical trajectories from the top. 

Traditionally, the equation for the dynamically corrected rate constant is derived by ap- 
plying a small perturbation to the equilibrium state and invoking the fluctuation-dissipation 
theorem and Onsager's relation [H[2tl25j. However, as I will show here, there is an alternative 
derivation that naturally evolves to a formula for transmission coefficient that is probably 
more efficient that the standard one |23p24]. 

There are several definitions for the rate constant k^B between two states A and B, such 
as the transition probability per unit time, the inverse mean residence time in state A, or the 
inverse mean first passage time towards state B [29]. However, all these different definitions 
become equivalent for truly exponential relaxation, which is the case whenever the stable 
states A and B are separated by large free energy barriers. If this is not the case, the rate 
constant becomes ill-defined. To start the derivation I will use the first definition, which can 
be expressed as follows: 

1 number of states A that transform into state B within dt 

kAB = lim (1) 

dt-s>o dt number of states A 

Let us denote x = (r, v) the phasepoint which includes the positions r and velocities v of all 
particles in the system. We define the reaction coordinate A(x) which can be any function of 
X, though in practice it will generally only depend on r. The RC function should describe the 
progress of the reaction, but there is a lot of flexibility in designing this RC function. 

We will assume that the collection of phasepoints {x|A(x) = 0} defines the transition state 
dividing surface that separates region A and B. For convenience, we will also assume that the 
RC will increase when going from A to B. Considering the phenomenological equation [H we 
can directly write down the reaction rate as 

, 1 Jdxoei-X{xo))6{X{x^t))p{xo) ^ 1 {6{-X{xo))e{X{x^t))) ... 

d™odt Jdxo9{-X{xo))p{xo) d™odt (^(-A(xo))) ^ ' 

where xq and x^t are phasepoints at times t = and t = dt. p{x) denotes the phasepoint 
density. For equilibrium statistics this is simply given by Boltzmann p[x) = exp{—f3E{x)) 
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where E the energy and /3 = l/ksT, T the temperature, and ks the Boltzmann constant. 9 
is the Heaviside-step function with 9{y) = if y < and 9{y) = 1 otherwise. The brackets 
(...) = / dx . . . p{x)/ J dx p{x) denote the ensemble average over the initial condition xq. Eq.[2] 
is basically the TST expression of the rate, but written in a somewhat unusual form. 

To transform this equation into the standard form, we can use A {xdt) = A (xq) + dtX{xQ) + 
C(dt^), where the dot denotes the time derivative. If we neglect the second order terms, we 
can write for an arbitrary function a{x): 

dxo e (-A (xo)) 9 (A {xAt)) a(xo) = j dx^O (-A (xq)) 9 (a (xq) + dtA (xq)) a(xo) (3) 



Clearly, 9{—X)9{\ + dtX) is only nonzero if A > and dt\ < A < 0. Instead of integrating 
over Xq, we will apply a coordinate transform such that we can integrate Eq. [3] over A, A and 
a remaining set coordinates Xg. Assume that J(xq, A, A) is the corresponding Jacobian of this 
transformation. We can then integrate out the (A, A) coordinates 



dxo ^ (— A (xq)) ^ (^A (xq) + dtA (xq)^ a(xo) = j dxg j dX j dA a(xo. A, A) J(xo, A, A) 
= / d4 f dA {(dtA)a(4, 0. A) J(4. 0, A) - 1 ,da)'^M2iiM^(i!kM) |^_^ + . . . } 

/roo 
dx'oj dAAa(x^,0,A)J(xo,0,A) + C'(dt2) (4) 

where we applied a Taylor expansion in terms of A in the second line. Clearly, as 

dxo A (xo) (5 (A (xo)) ^ (^A (xq)^ a(xo) = ^ dxgdAdA A5(A)6'(A)a(xo, A, A)J(xo, A, A) 

roo 

dxg / dA Aa(xo, 0, A) J(xo, 0, A) (5) 



we have proven that 

lim 1^ (-A (xo)) 9 (A (xdt)) = A (xo) 5 (A (xo)) 9 (x (xq)) (6) 
dt— >o at ^ / 

Using this expressing into Eq. [21 we obtain the standard form of the TST formula: 



{Xixo)5{X{xo))9lXixo)^. 
(^(-A(xo))) 

It is often very convenient to switch back to the other formalism, Eq. [21 as some relations 
follow more naturally from this expression, especially in path sampling simulations where dt 
can simply be taken as the MD timestep. The TST approach rewrites Eq. [7] into two factors 

_ (A(xo)c^(A(xo))g(A(xo))) (^(A(xo))) _ ^st , e'^^^^^ 

(HH^) (^(-A(xo)))-'' f.^dXe-m^) 
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where the free energy F is defined as F{X) = — ln(5(A — X{x)))/(3. Numerous techniques exist 
to calculate the free energy profile along the barrier region |27 tl281[30ll35j . The kinetic term 
^TST ^g^a,lly follows from a simple numerical or analytical integration. For instance, if the 
RC is a simple Cartesian coordinate of a target particle then R^^^"^ = 1/ y/2T:f3m where m is 
the mass of the particle. 

The TST expression neglects correlated fast recrossings and, therefore, overestimates the 
reaction rate. Recrossings can occur due to a diffusive motion on top of the barrier or by 
kinetic correlations when the kinetic energy of the RC is not dissipated. Another important 
source of recrossings is when the one-dimensional RC gives an incomplete description of the 
reaction kinetics [2j. To correct for recrossing we can apply the effective positive flux formalism 
which neglects the crossings that are not "effective". At each side of the barrier we define 
regions that are the stable regions A and B. These might be smaller than the regions that 
we associate to the product and reactant state. Entering ^ or i? implies that the system is 
committed to that side, i.e. it might leave region A oi B shortly thereafter, but the chance 
to rapidly recross the barrier is of the same order as an independent new event. An effective 
positive crossing is then defined as the first crossing on the trajectory that makes the transition 
from Ato B (See fig. [T]). This leads to the effective positive flux expression for the reaction 
rate: 
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Figure 1: Definition of an effective positive crossing on a very long MD trajectory. The EPF 
algorithm will ignore all crossing with the TS dividing surface except one (black dot). These 
are the first crossing points with the TS dividing surface for the parts of the MD trajectory 
that start at A and end at B (without revisiting A again). 



, _ y 1 {e (-A (xq)) e (A (xdi)) h%{xo)hi^ (xq)) 
- dt^odt (^(-A(xo))) 

(A (xo) S (A (xo)) e (a (xo)) h%{xo)hi^ (xq)) 

(^(-A(xo))) 



(9) 



where liJJ {xq) detects whether a backward/forward time trajectory crosses or enters a certain 
interface/region u before interface/region v. If this is true, the function is one. It is zero 
otherwise. In the second line we applied again equality El Naturally, the ensemble average 
(. . .) should now not only integrate over the phasepoint xq, but also sum over all possible 
trajectories backward and forward in time starting from xq. The ratio between the exact 
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expression (Eq. [9]) and the transition state expression (Eq. [7]) is the transmission coefficient: 



{X{xo)S{X{xo))e[x{xo)) 




{X{xo)6{X{xo))e 


(a(xo)) 


1) 



{x{xo)e(x{xo)) 


1 /iAo(^o)/lBA (3;o))a=0 


(A (xo) e (a (xo)) 


I)a=o 



(10) 

Here the subscript A = denotes an ensemble average on the TST dividing surface. Strictly 
speaking, the above expression is correct for any surface that separates the two stable states. 
However, the efficiency to calculate the above expression is significant better if k is maximized. 
Therefore, A = should be defined on the top of the free energy barrier. If we assume 
that A(x) = A(r) depends on configuration space only, the calculation requires to generate a 
representative set of configuration points on the TST surface. Then, we attribute to these 
points r a randomized set of velocities taken from a Maxwellian distribution and integrate the 
equations of motion backward and forward in time. However, as 6 (^X^ /i^q = if A < or 
when the backward trajectory recrosses the TST dividing surface before entering A, only a 
very few trajectories need to be fully integrated in both time-directions until reaching stable 
states. It is surprising that the effective positive flux counting strategy is not so common. 
To our knowledge only two slightly different expressions of a transmission coefficient based 
on the effective positive flux have been proposed in Refs. [361 137] . All other expressions in 
the literature do not avoid the counting of recrossings. In these algorithms, the final rate 
constant follows through cancellation of many negative and positive terms. For instance, the 
most popular formulation of the rate constant and transmission coefficient is the Bennett- 
Chandler (BC) expression that appears in many textbooks on molecular simulation [HE]. 

r (A(xo){(Afa))0(A(x.))) , 
= «*(-A(x„))> ^ 

~(.\ (A(xo)6'(A(xf)))A=o . . 

n(t) = — J-. r (11) 

(A(xo)e(A(xo))>A=o 



Here, the reaction rate and transmission coefficient are expressed as time-dependent functions. 
However, the actual rate constant and transmission coefficient, which should not depend on 
time, follow from a plateau value of these time-dependent functions: Uab = kAsit'), k = K{t') 
with Tmol < t' < Trxn- In Other words, k{t) and K{t) will generally show oscillatory be- 
havior at small t. However, after some molecular timescale t^oi, the system will basically 
enter either region A or B (See fig. [1]) after which we won't expect any recrossing until 
reaching the actual relaxation time Trxn >> Tmol- The equivalence between Eq. [11] and 
Eqs. [9l [lOl can be shown by invoking 9{X{xt')) = /i;gy^(xo), Eq. [6l and its mirror equiva- 
lent limM^(A(xo))0(-A(xdt))M = -X{xo)6{X{xo))e{-X{xo)): 



k = k{t' 



(A (xq) ,5 (A (xq)) e (A jxt'))) _ (A (xq) 5 (A (xp)) (xp)) 



(^(-A(xp))) (^(-A(xo))) 
(A (xp) 5 (A (xp)) e (a (xp)) 4^ (xp) + A (xp) 6 (A (xp)) 9 (-A (xp)) 4^ (xp)) 

(^(-A(xp))> 




-A (xq)) e (A (xdt)) 4a (^o) - (A (xq)) 9 (-A (xdt)) 4a (^^o)) \ ..21 
(^(-A(xo))) y ^ ^ 

We have now transferred the BC expression in an unitary ensemble average; each phasepoint 
Xq either returns 1, 0, or -1. Consider a very long MD trajectory with a timestep of dt (like 
the one in fig. [1]). It is clear that any detailed-balance simulation method should sample each 
phasepoint xq on this trajectory equally often. As such, an unreactive B ^ B trajectory will 
always have an equal number of phasepoints returning +1 as —1. The B ^ B trajectories 
are therefore effectively not counted due to this cancellation. The phasepoints on the A ^ A 
trajectory are always zero due to the /i^g^ characteristic function. Finally, any trajectory 
A ^ B always has one xq more that is +1 than -1. A more formal mathematical proof of the 
equivalence between Eq. [TT]and Eq. [U]can be found in Ref. |38] . 

Whenever, there are a significant number of recrosssings, the BC formalism has obvious 
disadvantages. In general, we note that any averaging method counting only zero and positive 
values will show a faster convergence than one that is based on cancellation of positive en 
negative terms. Moreover, in the effective flux formalism many trajectories will be assigned as 
unreactive after just a few MD steps, thus reducing the number of required force evaluations. 
Another important advantage of the EPF formalism is that is generates a set of trajectories 
that are unambiguous interpretable as reactive or unreactive, while the BC scheme generates 
only forward trajectories of which some actually belong to unreactive B ^ B trajectories. 
Instead of integrating the equations of motion until reaching stable states, one can also use a 
time-dependent expression for the EPF (39j similar to Eq. [TTJ 

There are several other formulations of the transmission coefficient (see [IDJ), but most 
of them rely on a cancellation between positive and negative flux terms. A comparative 
study of ion channel diffusion [51] showed that the algorithm based on effective positive flux 
expression was superior to the other transmission rate expressions. Moreover, it was as efficient 
as an optimized version of the more complicated method of Ruiz-Montero et al. [l2]. The 
implementation of the EPF scheme is as simple as algorithms that are based on the BC 
transmission coefficient. Therefore, the EPF implementation of the RF method should in 
principal be preferred above the standard implementations that require cancellation. 



3 Transition Path sampling 

In the previous section, I showed how the standard transmission coefficient calculations can 
be improved using the effective positive flux expression. However, this approach can not fully 
eliminate the main bottleneck of the RF methods. If k << 1 the number of trajectories that 
are required for sufficient statistics can be tremendous. In specific, if one is unable to find a 
proper RC, the overwhelming majority of trajectories that are released from the top of the 
barrier will be either A^ Aox B ^ B trajectories p]. In practice, it has been discovered that 
finding a good RC can be extremely difficult in high dimensional complex systems. Notable 
examples are chemical reactions in solution, where the reaction mechanism often depends 
on highly non- trivial solvent rearrangements [l3]. Also, computer simulations of nucleation 
processes use very complicated order parameters to distinguish between particles belonging to 
the liquid and solid phase. This makes it unfeasible to construct a single RC that accurately 
describes the exact place of cross-over transitions. As a result, hysteresis effects and low 
transmission coefficients are almost unavoidable. 
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This has been the main motivation of Chandler and collaborators [3H7] to devised a method 
that generates reactive trajectories without the need of a RC. This method, called transition 
path sampling (TPS), gathers a collection of trajectories connecting the reactant to the product 
stable region by employing a Monte Carlo (MC) procedure called shooting. 

Suppose X is a path {xq, Xdt,X2dt, • • • i Xndt} of n timeslices. The statistical weight given 
to this path equals 

P[x] = p{xo)pixo Xdt)p{xdt X2dt) ■ ■ ■ P{X(n-l)dt Xndt)H^) (13) 

where p{xo) is the usual phasepoint density and p{xjdt X(j+i)dt) is probability density that 
the MD integrator generates X(^j^i^dt starting from Xjdt- The characteristic function /i(x) 
equals 1 (otherwise 0) if a specific condition is fulfilled. For instance, one could imply that 
the trajectory x needs to start in state A and end in state B. 

By means of the shooting algorithm, TPS performs a random walk in path space to generate 
one trajectory after the other (See fig. [2]). The first step of this approach consists of a random 
selection of one of the timeslices of the old path, called the shooting point. This timeslice is 
modified by making random modifications in the velocities and/or positions. Then, there is 
usually an acceptance or rejection step based on the energy diff'erence between modified and 
unmodified shooting point. If accepted, the equations of motion are integrated forward and 
backward in time until a certain path length is obtained or until the condition function /i(x) 
can be assigned or 1. In the last case the trial move will be accepted. Any rejection along 
this scheme implies that the whole trial path will be rejected and the old path is counted again 
just like in standard Metropolis MC. Naturally, the random walk in path space should obey 




Figure 2: Illustration of the TPS shooting move using fiexible path length. From an existing 
path (black line) a random timeslice is selected. Positions and/or velocities of this point are 
slightly modified giving a new phasepoint (blue dot). From this point, the equations of motion 
are integrated forward and backward in time until the trajectory hits A or B. 



detailed balance 

Pgen[x(°) ^ x(°)] Pacc[x(°) ^ x^] _ P[xN] 
Pgen[x(-) ^ x(°)] Pacc[x(>^) ^ x(°)] P[x(°)] ' 



(14) 



where the superscripts (o) and (n) denote the old and new path respectively, and -Pgen[x — )■ x'] 
is the probability to generate path x' staring from x. Following the Metropolis-Hastings 
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scheme, the acceptance rule of the whole shooting move can be written as 



-face 



(o) 



X 



(n)l 



(n)^ 



mm 



' P[x(°)] Pgen[x(°) ^xH] 



(15) 



The generation probability is a product of different sub-probabilities. These are P^ei, to 
select the shooting point, Prarn for the random modification of this shooting point, and Ptrajj 
which is the probability to obtain x*^") by integrating the equations of motion backward and 
forward in time starting from the modified shooting point. 

-fgen[x y X ] = -Psel(^shoot|x) -fran (^^shoot ^ •^shoot) -^traj (x l^^shoot) (■'■^) 

If we generate paths of a fixed length n and if each timeslice has an equal probability to 
be selected then Pgei = 1/n. We come back to this point later on. In addition, TPS algo- 
rithms generally utilize a symmetric random modification of the shooting point: -Pran(3;shoot ~^ 
^shoot) = -Pran(2;shoot ~^ a^shoot)- Therefore, both Psei and Pran cancel out in Eq. [15i The ac- 
ceptance rule simplifies even further if we also assume that the dynamics obey the microscopic 
reversibility condition 

p{x)p{x -f y) = p{y)p{y x) (17) 

where x is the phasepoint x with reversed velocities: x = (r, —v). This relation is very general 
and valid for a broad class of dynamics applying to both equilibrium and non-equilibrium 
systems [7j . By applying Eq. [T7] several times on Eq. [13] we can show that 



P[x] = p{Xjdt)p{Xjdt X,^j__i)dt)p{X{j-l)dt X(^j_2)dt) ■ ■ ■ P{Xdt Xo) 
X piXjdt Xy+i^dt)pix(j+l)dt 3;(j+2)dt) • • • P{X(n-l)dt Xndt) 



(18) 



is true for any timeslice j. For time-reversible dynamics the backward integration is simply 
obtained by reversing the velocities and integration forward in time. Hence, the generation 
probability Ptraj(x|xshoot) depends on exactly the same transition probabilities p{x — )• x') and 
p{x —7- x'). This implies that all terms cancel out except the phasepoint density of the shooting 
point 



Pa, 



X 



(o) 



X 



Hi 



/i(x 



mm 



p[x^ ' 



shoot I 



P(2;Shoot^J 



(19) 



This is very convenient as this acceptance/rejection step can take place before the expensive 
trajectory generation takes place. Still, some (partly) completed trajectories will be rejected 
in the end due to the condition /i(x). However, as the new trajectory was generated from a 
small modification of an existing trajectory with /i(x) = 1 the chances are relatively high that 
the condition will be satisfied for the trial trajectory as well. 

The sampling of trajectories under a given condition h might benefit from using a path 
ensemble that has a non-fixed length. Using a fixed path length to sample all possible tra- 
jectories between A and B is expensive as this length needs to be adapted to the longest 
pathway connecting these states. Many trajectories will reach A to i? in a much shorter time 
and will, therefore, consist of unnecessary parts that are not relevant for the actual barrier 
crossing event. In addition, if trajectories have significant parts outside the barrier region, 
the shooting move becomes inefficient as many shooting points will lie inside the reactant or 
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product well. Shooting from these points gives a very low probability to connect both states. 
Using flexible path lengths was first introduced in Ref. [8] within the context of the TIS rate 
evaluation. However, also for the generation of reactive trajectories, the flexible path ensem- 
ble is very useful and allows to generate paths that start and end just at the boundaries of 
A and B (see Fig. [2]). The only difference with the previous example is that Pgci = 1/n is 
not cancelled as the trajectory length can be different. Therefore, if the shooting procedure 
selects the timeslices by an equal probability, the acceptance rule becomes 



-face [■ 



X 



(o) 
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(n)l 



/i(x 



mm 



in) . 
shoot/ 



n 



(o) 



(o) \ 

shoot / 



(20) 



with n^°\n^^^ the length of the old and new path. This expression is not so convenient as 
a rejection can only be made whenever the whole path is completed. Hence, the integration 
needs to be carried out even if p{x^2!oot) /'(^shoot) implying an almost certain rejection. 
We can, however, separate the acceptance into two steps by writing 



-face [ 
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(o) 
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(n)l 
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(n)^ 



mm 



/^(^shoot) 
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X mm 
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(n) 



(21) 



This acceptance rule obeys detailed balance as well and allows to reject a modification of the 
shooting move that gives a too high energy. Still, even if the first step is accepted, the final 
trajectory might be rejected whenever it becomes too long compared to the previous path. We 
can improve the efficiency even further using following trick [8j. Instead of taking a random 
number a G [0 : 1] after finishing our trajectory and then accept if a < n^"^/n^^\ we will 
actually draw this random number before starting the integration of motions. As we now 
know that we will have to reject our trajectory whenever a < n^°^ /n^") , we can simply define 
a maximum allowed path length of this trial move in advance 



n 



int [n(°Va] 



(22) 



This allows to directly stop our trial move whenever it exceeds this maximum path length. 

The original TPS method also provided an algorithm to determine the reaction rate of the 
process. This approach has been improved by the TIS |8j and RETIS [^llOj algorithms. Like 
RF, the TPS rate evaluation does require a RC (I will not make the distinction between order- 
parameter or RC). However, one can show that, compared to the RF method, the efficiency 
of TPS, TIS, and RETIS, is less sensitive to an improper choice of the RC [39]. 

The original TPS rate evaluation is based on following correlation function 



C{t) 



{hAixo)hB{xt)) 
{hA{xo)) 



(23) 



where hA/Bi^) = 1 if G A/B and otherwise. Just like Eqs. [TT| C{t) will initially show 
some oscillations. However, if there is a separation of timescales, this correlation function 
grows linearly in time, C{t) ~ kABt, for times Tmol < ^ < Trxn- Hence, 



d (hAixo)hB{xt)) 

kAB{t) = -C{t) = ^ ,, , ^ kAB = kAB{t') for 

dt {hA[xQ)) 



Tmol <t < Trxn (24) 
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The correlation function C{t) is calculated in the TPS scheme using the shooting algorithm in 
combination with umbrella sampling. First, the fixed path length t' is fixed to a value where 
C{t) should give a plateau. Then a series of path sampling simulations will be performed 
in which the final region B is slowly shrunk in successive steps from the entire phase space 
to the final stable state B [7j. For each step numerous trajectories are generated with that 
condition that the path should start in A and end in the extended region B at time t'. 
The distribution of the path's end-point will be binned into histograms that will be matched 
just like ordinary umbrella sampling. Once the fully matched histogram is obtained, C{t) is 
obtained by integration of this histogram over the actual region B. 

The approach is rather time-consuming because it can take a relatively long time Tmoi 
before C{t) reaches a plateau (longer than in a transmission coefficient calculation [7]). In 
Ref. [5] an improvement of this approach was presented in which the umbrella sampling series 
could be performed with paths shorter than t^oI- The results were then corrected by a factor 
that is obtained from a single path sampling simulation using the longer paths. Unfortunately, 
the relative error in this correction factor is large if the path length is reduced too much, so 
that the gain in CPU efficiency remains limited [8]. Moreover, inspection of Eqs. ()23p and 
(j24|) shows that a necessary cancellation of positive and negative terms can slow down the 
convergence of the MC sampling procedure. 



4 Transition Interface Sampling 

TIS is a more efficient way to calculate reaction constant than the method discussed above. 
The TIS methodology is also the basis of several other algorithms [9lll H[T3ll44p 45j of which the 
PPTIS, RETIS, and FFS methods will be discussed in forth-coming sections. The TIS rate 
equation is related to both the EPF expression, Eq. [9l and to the correlation function used 
in TPS, Eq. [23l albeit using different kind of characteristic functions. Instead of using the 
characteristic functions of the stable states A and B, we will redefine the correlation function 
using overall states A and B. These states do not only depend on the position at the time of 
consideration but also on its past behavior. Overall state A covers all phase space points lying 
inside stable region A, which constitutes the largest part, but also all phase space points that 
visit A, before reaching B when the equations of motion are integrated backward in time. In 
other words, all phasepoints that were more recently in A rather than in B. Similarly, state 
B comprises stable state B and all phase points, coming directly from this state in the past, 
i.e. with- out having been in A. The corresponding correlation function is 

C{t) = (Mxo)M^O) (25) 

where /i^(xo) = h\Q{xo) and Hq^xq) = /i^^(xo). Contrary to Eq.[231 this correlation function 
has no oscillatory behavior during a molecular timescale Tmoi- On the contrary, it exhibits a 
linear regime ~ k^st for < t < r^xn- The system will only transfer from overall state A to 
overall state B when it enters region B for the first time since it left region A. If it leaves 
state B shortly thereafter, it will remain in B. Therefore hQ{xt) and hji,{xt) do not show the 
fast fluctuations that are found for hB{xt) and hA{xt)- As Eq. [25]is linear from the start, we 
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can simply take the time derivative at t = 0, which gives 



kAB = - — n — 7 — — = hm ^ 



(/i^(xo)) dt->o dt (hAixo)) 

h\j,{xo)X (xo) 6{X{xo) - XB)e (X {xo)] 



(26) 



The resulting expression is basically the EPF expression (Eq. [9]) through the interface A^- 
At a first sight it might seem that generating long trajectories is no longer needed. As we 
only need dC{t)/dt at t = 0, the minimum time range over which we need to calculate C{t) 
is [0 : dt] instead of [0 : Tmoi]- Unfortunately, unlike hAixo) and hsixo), the determination 
of hj^ixo) and hig(xo) can not be done instantaneously. For this we still need to integrate 
the equations of motion. However, for most xq, hA{xo) /hB{xo) can be assigned 1 or using 
a much shorter backward trajectory than Tmol- For stochastic dynamics h^ixQ) /hQ{xQ) can 
be, strictly speaking, a fractional number. However, there is generally no need to know this 
fractional number for a specific phasepoint, except for committor analysis |14ll46fH9] . Hence, 
TIS algorithms will generally compute /iyi(xo)//i_B(xo) for one specific path to which xq belongs. 
Conceptually, it is therefore more accurate to speak of a MC sampling in pathspace rather 
than phasespace. The TIS correlation function has an additional advantage that the reaction 
rate is somewhat better defined if the separation of timescales Tmol ^ '^rxn is not sufficiently 
obeyed. The fact that the derivative of C{t) is taken at t = makes corrections like the one 
suggested in Ref. [SD] unnecessary. 

The TIS algorithm expresses the rate equation, Eq. [26l as a product of different terms. 
Each term has a much higher value than the final rate and is, therefore, much easier to 
compute. To introduce the TIS and PPTIS expression, that I will discuss in the next section, it 
is convenient to introduce following crossing probabilities that depend on four non-intersecting 
interfaces {x|A(x) = Aj},{x|A(a;) = Aj},{x|A(x) = Xk}-,{x\X{x) = A;} 

r {h%{x,)e{X, - A(xo))g(A(xd,) - A,)/.{,(xo)) 



dt^o (/i^^.(xo)0(A,- - A(xo))0(A(xdO - A,)) 

(/i^^.(xo)A(xo)<5(A(xo)-A,)e(A(xo))4(xo)> ^ ^ 

= ^ — 7 — - — r for A,- > Xi (27) 

{h%{xo)X (xo) 5 (A (xo) - A,) e (a (xo)) ) 

For Xj < Xi we simply need to replace 9{Xj — X{xo))9{X{xdt) — Xj) by 9{X{xq) — Xj)9{Xj — X^x^t)) 
in the first line or A (xo) 6 (A (xo) — Xj) 9 (^X (xo)^ by —A (xo) 6 (A (xo) — Xj) 9 ^— A (xo)^ in the 
second line of the above definition. Eq. [221 defines a conditional crossing probability. It is the 
probability that the system will cross interface A^ before A; under a twofold condition. These 
conditions are that the system should cross interface Xj right now at time t = 0, while Xi was 
more recently crossed than Xj in the past. 

Using these crossing probabilities, one can prove that Eq. [26] is equivalent to the product 
of the initial fiux times the overall crossing probability [8j 



(A(xo)5(A(xo)-Ao)e A(xo) ) 

kAB = , ^ ^ X P(X-) ^ fAVA{Xn\Xo) (28) 
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where Aq and are the boundaries of the stable states A and B. /a is just the flux out of state 
A that can be computed with standard MD as the boundary of A is set at the left side of the 
barrier region. The minus in 0~ is to denote an interface Aq — e which is put there to indicate 
the direction of the crossing at t = 0. The overall crossing probability Pyi(A„|Ao) = -P(olo-) is 
the probability that once Aq is crossed, A„ will be crossed before a recrossing with Aq occurs. 
This probability is very small, but it can be calculated by defining n — 1 non-intersecting 
interfaces in between Aq and A^ and express the overall crossing probability as the following 
product [8] 

n-l 

7^a(A„|Ao) = n^A(A.+i|Ai) (29) 

The factorization of ^^(AnlAo) into probabilities VAiK+i\K) that are much higher than the 
overall crossing probability, is the basis of the importance sampling approach. It is important 
to note that VAi^i+iW) are in fact complicated history dependent conditional probabilities. 
If we consider all possible pathways that start at Xa and end by either crossing Xa or As, 
while have at least one crossing with Aj in between, the fraction that crosses Aj+i as well 
equals 7^A(Aj+i| Aj). This basically reduces the problem to a correct sampling of trajectories 
that should obey the Aj crossing condition (See fig. E]). From now on we will call this the 
[i^] path ensemble. In TIS, this is done via the shooting algorithm for flexible path lengths 
as is discussed in previous section (For a full flowchart diagram of the TIS algorithm see 
Ref. [51]). The number of interfaces and their separation should be set to maximize efficiency. 




Figure 3: The TIS path ensemble [i"*"] is required to calculate the conditional crossing prob- 
ability 7^^(Aj+i|Ai). For this purpose we apply the shooting move to generate all possible 
trajectories starting at Aq and ending at Aq or A^ with at least one crossing with Aj. Suppose 
the algorithm starts with the middle path that already fulfills these requirements. A shooting 
point is randomly selected and modified (black dots). However, the trajectory that starts out 
from this point (top trajectory) fails to cross A, and is therefore rejected and the old path is 
counted again. A new shooting (blue dots) generates a valid trajectory that not only crosses Aj 
but Aj+i as well. This trajectory is called "successful". The fraction of successful trajectories 
in this ensemble determine P^(Aj+i|Aj). 

In Ref. [39y4U) it was found that the optimal interface separation is obtained when one out of 
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five trajectories reach the next interface. In addition, one can define a set of sub-interfaces of 
arbitrary separation in order to construct the crossing probability as a continuous function. 
This strictly decreasing function could be viewed as the dynamical analogue of the free energy 
profile F(A). 

There are some small differences how to treat the path ensemble [i"*"] regarding the end- 
point of the path. In the first TIS algorithms, the trajectory could reach upto Aj+i where 
the trajectory was stopped and assigned successful. More recent simulations continue the 
trajectory until reaching the stable states A oi B each time. The additional cost is very 
limited as about 80% is not reaching Aj+i and need to be followed until reaching A anyway. 
The choice to continue the trajectory even after Aj+i has the advantage that one can start the 
[i^] path ensemble without the need to fix a value for Aj+i beforehand. After some simulation 
cycles the Aj+i can be set to have the optimal 20% success-rate after which one can start 
the [{i + 1)^] path ensemble. In addition, the new approach makes it much more easy to use 
replica exchange which we will discuss in Sec. [71 

The simplicity of Eq. [2U]is deceptive and could be mistaken as a Markovian approximation. 
The reason that the equation is still exact lies in the fact that the crossing probabilities are 
history dependent and by the fact that it only considers first crossing events. We can argue the 
exactness of the equations also in another way. Suppose we want to calculate the probability 
to go from Ao to Ai to A2 . • • to A„ in successive jumps. This probability can be expressed as 

P(Ao^Ai^A2^...^A„) = P(Ao)xP(Ao^Ai|Ao)xP(Ai^A2|Ao^Ai) 

X P(A2 A3IA0 Ai A2) X . . . 
X P(A„_i A„|Ao ^ Ai ^ A2 ^ . . . A„_2 A„_i)(30) 

This is an exact non-Markovian expression for this specific crossing sequence that looks similar 
to Eq. [29j However, it doesnot say anything about the many different trajectories that could 
connect Aq with A„. For instance, we should also take into account the sequence Aq — )• Ai — >■ 
A2 — ^ Ai — )• A2 — )• A3 — )• . . . A„_i — )• A„. Therefore, it might seem that the right expression 
should look much more complicated than Eq. [29l The trick, however, is that this last sequence 
can not occur if we only consider first crossing events. When we move back to Ai in the third 
step, this move will simply not considered as it is a second visit since leaving Aq. Hence, the 
successive sequence Aq — A2 — >■ A3 — )■ . . . A^-i Xn is the only possible sequence of first 
crossing events that brings you from Aq to A„. 



5 Partial Path Sampling 

The PPTIS is a variation of the TIS algorithm that was devised to treat diffusive barrier 
crossings Despite the existence of a fine separation of timescale, i.e. the time to cross 
the barrier is still negligible compared to the time spend in the reactant well, the path length 
can become too long for an effective computation of the reaction rate. This is the case if 
the barriers are sufficiently high to ensure exponential relaxation, but not very sharp so that 
system can move backward and forward on the barrier before it eventually drops off. The 
PPTIS equation depends on the same rate equation as TIS 

kAB = fAPQl) (31) 
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The only difference witli Eq. [28] is that we now consider the condition |q) instead of |g-), but 
this is just a technical detail. If we take Ai = Ao + e the equations become equivalent. However, 
Ai can be any value that is somewhat larger than Aq by redefining Ja as the effective flux 
through Ai. This implies that we should count the positive crossings with Ai whenever the 
system leaves the stable state A : {x\X{x) < Aq}- However, the next positive crossing should 
only be counted if the system has revisited A again. The PPTIS approach tries to avoid 
the generation of very long trajectories using a soft Markovian approximation. The PPTIS 
scheme assumes that for a well positioned set of interfaces the system will lose its memory 
over a distance that is similar to the interface separation. This implies for any m > 1 

Pik\'j±m) ^ -f(Llj±l) (32) 

The PPTIS algorithm consist again of a series of path sampling simulations. Each PPTIS 
simulation samples a certain path ensemble in which trajectories are confined within two next- 
nearest interfaces. For instance, the [i^] path ensemble will consist of all possible trajectories 
starting and ending at either Aj+i or Aj_i having at least one crossing with the middle interface 
Aj. From these simulations the following short-distance crossing can be obtained 

pt^p(twi-i), pr^p&-i), pr^p(xM+i), A^pQ^\i+i) (33) 

with pf + Pi" = pf + pI = 1- For instance, pf is determined by dividing the number of 
trajectories in the [i^] ensemble that start at Ai_i and end at Aj+i divided by all trajectories 
that start at Aj_i. 

Once these short distance crossing probabilities are obtained with sufficient accuracy, the 
overall crossing probability can be obtained. One way to do this is to use these probabilities as 
input for a kinetic MC simulation [52]. However, this is not needed for a one-dimensional RC 
which allows an elegant analytical treatment [11]. Naturally, -P(olo) — P^^ but the calculation 
of -P(olo) requires already to sum up the trajectories 0— t-I— )-2— )>3, 0— )-l— )-2— )-l— )'2— >-3, 
etc. However, as shown in Ref. jllj . one can derive following recursive relations to make this 
infinite summation of all trajectories (include the ones of infinite length!). These PPTIS 
recursive relations are the following 

^± p(m\l\ \m—l\ 
T3(m+l\l\ "m VP 10/ P/'O 1"^ ^ rm \m\m J ('iA\ 

Pm^Pm^\ra\m ) Pm^Pm^\m\m ) 

or, by defining the long-distance crossing probabilities = P(™|q),P^ = -P(mlm^^) 
kAB = JaPu 

^nl+i = ^fp- ^ P^+i= /f?p- . form>l, Pt = P{ = l (35) 

Hence, starting from the initial conditions for {P^,P^) = (1,1), one can successively solve 
(^2^ ) P2 ) ■> {P^ ' -^3" ) 5 • • • ) {Pn 1 Pn ) '^i^ Eq. [35l It is important to note that P~ is not exactly 
the same as P^ in the reverse direction. Only for j = n these two probabilities can be viewed 
as mirror images. 

Here, I will derive an alternative recursive relation that does not require the auxiliary 
reverse probabilities P^ . The derivation is similar to the one presented in the supplemental 
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information of Ref. |53j which treats the simpler hopping process. To achieve this, we will 
bring P{m\m~^) in front of Eq. 



Pi 



im— I'l 
m\m J 



or by incrementing m 



PC 



m+1 Im+l; 



(36) 



(37) 



Moreover, we can write for ^'(^'''^lo; 



pfm+l\l\ _ p/m|l\ p/m+lim ^ 
\0) — loi-'^lo \m-l) 



-P((rio) (l ~ -f(m+llm-l)) 
Pio^lo) (l ~ -f (m+llm+l)Pm/Pm 



Then, we substitute Eq. [37] in Eq. [38] and bring -P(o^'''^|o) in front which yields 



p/^n+2|l^ 



p/m+i |l% 



(38) 



(39) 



or 



kAB 



fAP;l 



PmPm+ 1 Pm Pm+ 1 



m+2 



+ PmPm+l] P" 



PrnPm- 



In the case of a full Markovian assumption, pf = p\ and = p^^, we reobtain the simpler 
expression of Ref. |53J. For some this new expression might seem more esthetic as the set 
of two equations has now been transferred into a single recursive equation without relying 



Pt 



1, 



Pt 



(40) 



on the auxiliary probability P^. The new function has its utility 



but is numerically 



somewhat problematic as it can sometimes produce zeros in both nominator and denominator 
that cancel, but is not practical for numerical calculations. 

The positioning of interfaces is crucial in PPTIS. On one hand, one would like to put them 
close together to improve efficiency. On the other hand, putting them to close will introduce 
systematic errors due to a decrease of history dependence of the hopping probabilities, which 
invalidates Eq. [32] A way to measure whether the interfaces are sufficiently far is the calcu- 
lation of a memory-loss function However, the memory- loss function can only provide 
a necessary but not necessarily sufficient condition for this separation. In fig. [3] we give two 
examples of well and badly placed interfaces. 

The milestoning method [.12j is very similar to PPTIS. There are basically two important 
differences. Milestoning assumes full memory loss once the system hits an interface. In our 
notation this could be rephrased as P(^|^-^^) ~ Pikl])- This is a stronger approximation 



than Eq. [32] The approximation of milestoning becomes exact, if the interfaces coincide 
with the iso-committor functions [33], but these are difficult determine. On the other hand, 
milestoning is more precise in the construction of the time-evolution of the system by making 
the crossing probabilities time-dependent. This is important if there is not a clear separation 
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of timescales and also allows to calculate other dynamical properties like diffusion. Hence, 
instead of ,pT ^pf :Pi with pf + Pi = + p\ = ^, milestoning calculates for each interface 
the time-dependent probability densities pt{t),p^{t) with J^\pf{t) +p^{t)]dt = 1. The 
strengths of PPTIS and milestoning do not exclude each other and could be unified into a 
single method as was suggested in Ref. [55]. A realization of such a method was recently 
published 



Au Au Au Au Au Au Au Au 

1 2 3 1 2 3 4 




Figure 4: Examples of well and badly positioned interfaces with respect to the memory-loss 
assumption Eq. [32j The top situation requires a good description of the kinetic correlations 
whenever the top of the barrier has many small local wells. At the left, the interfaces are 
correctly placed. Once the system crosses Ai it will gain a lot of kinetic energy when arriving 
at A2- Henceforth, the chances are high that the system won't get trapped and directly moves 
upward to cross A3. The PPTIS simulation for this interface configuration will show that 
^> ^ as it should. At the right-hand side we have put an additional interface inside the local 
well. The [3^] path ensemble will consists of trajectories having a much lower kinetic energy 
than the [2^] ensemble of the left-hand side. Henceforth, the right-hand side will overestimate 
the probability to get trapped. The bottom picture shows impermeable wall (thick black line) 
with a small hole. The left hand side shows a correct positioning of interfaces. The pathways 
that are generated from Aq to A2 all have to move through the small hole. Conversely, the 
right-hand side shows a bad overlap between the [1^] and [2^] path ensembles which might 
give the impression that trajectories can tunnel through the wall. 



6 Forward Flux Sampling 

EES was originally developed for the special case of biochemical networks that do not obey 
equilibrium statistics nor time-reversibility [13j. However, its advantageous implementation 
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and apparent efficiency has gained this method a fast increasing popularity for equihbrium 
systems as well. FFS is based on the same theoretical TIS rate equation [26j However, the 
fundamental difference is the sampling move. While the principal sampling move in TIS is the 
shooting move, FFS is based on a non- Metropolis MC scheme called splitting |57tl58j. This 
approach requires stochastic dynamics, although it has been suggested that FFS is able to 
treat deterministic dynamics utilizing the Lyapunov instability [59''60] using small 'invisible' 
stochastic noises [61]. Like TIS, FFS consists of a straightforward MD simulation, from which 
the escape flux is obtained, followed by a series of path sampling simulations. However, 
besides giving the flux value, the MD simulation also provides the starting conditions for the 
path sampling simulations. Each time that the first interface Aq is crossed in the positive 
direction, this phasepoint just after the interface is stored on the hard-disk. In the first path 
simulation, performing the [0"*"] path ensemble, these points are used as starting points for the 
trajectories that are continued until reaching Ai or returning to Aq- Naturally, stochasticity 
is of eminent importance, otherwise all these trajectories would just reproduce parts of MD 
simulation. The endpoints of the trajectories that successfully reach Ai are stored again and 
serve as initial points for the [1^] ensemble. The path ensembles are executed one after the 
other by the same procedure until reaching state B (See fig. E]). 




Figure 5: Illustration of FFS algorithm, a) shows the initial MD simulation that is needed to 
calculate the flux /a- The positive crossing points (black dots) are stored, b) Starting from 
the stored MD crossing points a number of trajectories are released. The endpoints of the 
successful trajectories (that reach Ai) are stored again and used for the next path ensemble 
simulation c). Finally, the reactant state will be reached d). 

There are advantages and disadvantages compared to the TIS algorithm. The most impor- 
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tant advantage of FFS is that it allows to treat non-equilibrium systems as it does not require 
any knowledge about the phasepoint density. TIS employs the shooting move that requires 
to know p{x) for the acceptance, Eq. [T9j In addition, FFS does not require any integration 
of motion backward in time. Therefore, time-reversibility is not required. Moreover, unlike 
PPTIS, FFS is, in principle, equally exact as TIS. However, if one has to chose between TIS 
and FFS for equilibrium dynamics, one has to consider following points. FFS will generally 
create more trajectories for the same number of MD steps as it recycles previously generated 
trajectories. Moreover, there are no rejections like there are in TIS and any other Metropolis 
based MC scheme. In practice, the reduction in MD steps will be limited to a certain factor 
(ps 2,3) as the unsuccessful trajectories, which is the largest part, have to be followed until 
reaching state A. On the other hand, the FFS trajectories will be much more correlated 
than the TIS trajectories. This implies that FFS needs much more trajectories than TIS 
to obtain the same accuracy. One reason for this, is that FFS generates several trajectories 
having the same starting point. Absence of stochasticity will result that all these trajectories 
basically coincide. Fully Brownian motion does not exclude correlation effects either as the 
successful trajectories starting from the same point will hit the next interface in a confined 
region. The size of this region is determined by the diffusion orthogonal to the RC and the 
time it takes to go from interface to the other. Besides correlations within a certain path 
ensemble, the FFS method also introduces correlations between the different ensembles. This 
is a crucial difference with TIS where the MD simulation and all path simulations are inde- 
pendent. One of its consequences is that the FFS is more sensitive to the RC than TIS or 
even RF [39]. An efficiency analysis of FFS [62] ignores this correlation effect. This can be 
a rather crude approximation that is probably only valid when interfaces are approximately 
equal to the isocommittor surfaces. Suppose A"*" denotes a coordinate orthogonal to the RC. 
Let VA{^n\^o'i ^'^) be the overall crossing probability from Aq to A„, starting from a point A"*" 
on the first interface. Then, the full overall crossing probability is given by 



where p(A |Ao) is the probability density of A on interface Aq. FFS will suffer considerably 
when the distributions p(A-'-|Ao) and 'PA(An|Ao; A"*-) are not overlapping. In that case, FFS 
will miss important crossing points that are significant for the rate evaluation. Some studies 
have shown that FFS can significantly underestimate reaction rates |63y 64) in practical cases. 
Sampling artefacts like this, are also not yet fully excluded as possible explanation for some 
surprising results on non-equilibrium nucleation [65j . 

This issue will be most sensitive to the MD and the first interface ensembles on which all the 
further results will depend. If Ai, A2, . . . , A„_i are isocommittor surfaces then VA{^n\K', A"*") 
is a constant as function of A^ which eliminates the problem. This is the reason that Borrero 
et al. devised a FFS scheme in which the interfaces are repositioned on-the-fly in order to 
obtain a proper RC [66j. 

TIS has the advantage that it can relax the history of the path via the backward integration. 
Therefore, the distribution density p{X-^\Xo) can change when considering the different path 
ensembles. TIS can give correct results even if the sampled distribution of starting points A"*" 
of the final [n — trajectories do not overlap with the initial MD crossing points [67j . 




(41) 
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7 Replica-Exchange TIS 



In Ref. [9], I showed that a special type of repHca exchange |68ll69j can significantly improve 
the TIS efficiency (See also Ref [lOj for some extensions of this approach). A crucial difference 
with standard RE, which has also been applied to TPS [70j, is that the RETIS method does 
not require additional simulations at elevated temperatures. Instead, swaps are attempted 
between the different TIS path ensembles. For this purpose, RETIS has replaced the initial 
MD simulation by another path ensemble, called [0~], that consists of all path that start at 
-^0 = ^A, then go in the opposite direction away from the barrier inside state A, and finally 
end at Aq again. The flux is then obtained from the average path length of the [0~] and [C*"] 
ensembles as follows [9]. 



where {tp^^y, (^path) average path lengths in the [0^] and [0^] path ensembles respec- 

tively. 

As the dynamical process is now fully described by path simulations with different interface- 
crossing conditions, the exchange of trajectories between them becomes extremely efficient, 
especially if the process possesses multiple reaction channels [67]. The methodology avoids 
the need of doing additional simulations at elevated temperatures and even gives paths for 
free as for most swapping moves whole trajectories are being swapped. Only when a swapping 
between the [0~] and [O"*"] ensembles are attempted, two phase points are interchanged. From 
the last point of the [0~] trajectory a new path in the [0"^] is generated. Reversely, the first 
point of the old [0"*"] path will serve to generate a new path in the [0~] ensemble by integrating 
the equations of motion backward in time (See Fig. [6]). 

The RETIS algorithm is then as follows. At each step it is decided by an equal probability 
whether a series of shooting or swapping moves will be performed. In the first case, all 
simulations will be updated sequentially by one shooting move. In the second case, again an 
equal probability will decide whether the swaps [0~] -H- [0^],[1+] -H- [2+],... or the swaps 
[1+] o [2+],[3"'"] ^ [4:~^],... are performed. Each time that [0~] and [{n - 1)+] do not 
participate in the swapping move they are left unchanged. Also when the swapping move 
does not yield valid paths for both ensembles, the move is rejected for these two simulations 
and the old paths are counted again. Note that the swapping moves do not require any force 
calculations except for the [0~] -H- [0+] swap. 

Like EES, the path ensembles in RETIS are not fully uncorrelated. However, their de- 
pendence fundamentally different. In EES, the path ensemble is fully determined by its 
predecessors, the MD simulation and the path simulations [j"*"] with j < i. Conversely, the 
separate RETIS simulations generate a large part of their trajectories independently. The ex- 
change between the ensembles is, therefore, an additional help instead of a strict dependence 
as it is for FES. Moreover, the benefit of the exchange works in both directions and is mutual 
for all ensembles, i. e. the [i^] path ensemble can improve the sampling in both [{i + 1)^] and 
[{i — 1)^] via the swapping moves -H- [{i — 1)+] and [i^] -H- [{i + 1)^], and will improve 
itself due to the same moves. 
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Figure 6: Parallel path swapping move in RETIS. The picture illustrates a RETIS simula- 
tion using four interfaces, a) shows the initial "superstate" that contains one trajectory per 
ensemble [0~], [0+], and [2+]. b) shows the superstate after the [0+] o [1+] swap. The 
original [0"^] trajectory crossed Ai and is therefore a valid path in the [1"^] ensemble. The 
swapping move is, therefore, accepted, c) shows the trial superstate that is obtained after the 
simultaneous swaps [0~] -H- [O"*"] and [l"*"] -H- [2"'"]. The first swap requires the integration of 
motion forward or backward in time starting from the last of first timeslice of the swapped 
trajectories. This swap will always generate acceptable trajectories for the [0^], [0^] ensem- 
bles. The other swap [1+] o [2+] is rejected because the old [1+] path does not cross A2. e) 
gives the final situation after the whole move. 



8 Numerical Example 

We will apply the different methods on a simple one-dimensional test system using Langevin 
dynamics with finite friction. The Langevin dynamics was chosen because its inhibits stochas- 
ticity which is required for FFS. Hitherto, most studies on FFS have applied Brownian dy- 
namics. As the dynamics we are considering are not overdamped, the dimensionality of the 
system becomes effectively two-dimensional. We could consider the velocity as an orthogonal 
coordinate (like in Eq. HT]) . This makes the choice of a proper RC not such a triviality 
as one would think at first sight. However, we will simply take the RC to be configurational 
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dependent, which is the standard approach. The system that we will consider consist of a 
single one-dimensional particle inside a double well potential 

V{r) = kir'^ - fcar^ (43) 

with A;4 = 1 and k2 = 2. The corresponding potential has a maximum at r = and two 
minima at r = ±1. We use reduced units where the mass and the Boltzmann constant are 
set to unity, ks = m = 1. The system is coupled to a Langevin thermostat with friction 
coefficient 7 = 0.3 and temperature T = 0.07. The equations of motion are integrated using 
MD timestep of dt = 0.002. The RF method was applied using the EPF formalism for the 
transmission coefficient calculation. For this purpose 100,000 trajectories were released from 
the TST dividing surface r = 0. The free energy term dA e~^^^'^'^ was obtained by a 
simple numerical integration. The RF method is by far the most efficient method for this 
system because there is basically no error in the free energy calculation and the transmission 
coefficient is close to unity. The RF results will therefore be the reference for the other 
methods. The escape flux Ja for TIS, FFS and PPTIS was determined using a MD simulation 
of 10,000,000 timesteps. The same MD result was used for these three methods. We performed 
an additional MD simulation using less timesteps, 4,000,000, for a second FFS calculation in 
order to see how this effects the final FFS result. I defined eight interfaces Aq = —0.9, Ai = 
-0.8, A2 = -0.7, A3 = -0.6, A4 = -0.5, A5 = -0.4, As = -0.3, and A7 = 1.0. For each path 
ensemble 20,000 trajectories were generated. For TIS and PPTIS, 50% of the MC moves 
were shooting moves. I applied the aimless shooting [Tlj approach in which the velocities at 
the shooting point are completely regenerated from Maxwellian distribution. However, unlike 
Ref. |71) . shooting points were picked with an equal probability for all timeslices along the 
path without considering the previous shooting point [5]. The other 50% were time-reversal 
moves. Time-reversal moves simply change the order of the timeslices of the old path while 
reversing the velocities. Time-reversal can sometimes increase the ergodic sampling and is 
basically cost-free as it doesnot require any force calculations. However, as aimless shooting 
is also able to reverse velocities in a single step, the time- reversal move could actually have 
been omitted for this case. In the RETIS algorithm there was at each step a 25% probability 
to perform a shooting move, another 25% probability to do a time-reversal move, and a 50% 
probability to do a replica exchange move. The FFS simulations consist of a single move 
which is the forward integration of the equations of motion. The Langevin thermostat served 
for the necessarily stochasticity. The results are shown in table [H The RF method gives the 
most accurate results as expected. The value for k can be compared to Kramer's expression 
K ^ {l/ujb){-j/2 + ^72/4+7^) = 0.9 with ujb = y/h/m = \/2. If we compare the TIS, 
PPTIS, and RETIS simulations we see that they are all close (within a factor 2) to the RF 
result. The TIS result is somewhat too high which is probably due to a single path ensemble 
calculation that wasn't fully converged after 20,000 steps. The RETIS results are clearly 
better despite the apparent errors that are somewhat smaller for TIS. The RETIS result is 
much closer to the RF reference. Moreover, it uses only halve the number of shooting moves 
compared to TIS, which is the most expensive move for realistic systems as it requires a large 
number of force evaluations. Also, the construction of the overall crossing probabilities in 
Fig. El shows a much better matching between the different ensembles in the RETIS method. 
The PPTIS result is also very close to the reference value. The PPTIS approximation, Eq. [321 
becomes not only exact for very diffusive systems, but is also exact for steeply increasing 
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Table 1: Results of the rate evaluations using RF, TIS, PPTIS, RETIS, and FFS. Final errors 
were obtained by block averaging and error-propagation rules. The errors of RETIS and FFS 
are given a star as these errors should not be considered exact due to the neglect of covariant 
terms which arise due the correlations between path ensembles and initial MD simulation. 
The FFS was repeated using a shorter (4,000,000 instead of 10,000,000 timesteps) initial MD 
run. 
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Figure 7: Overall crossing probabilities using TIS, RETIS, and FFS 



barriers as all trajectories from Aj to Aj+i come directly from Aq in the past. FFS on the other 
hand, that is in principal exact unlike PPTIS, gives an unacceptable value that is about a 
factor 20 too low. Still, if we calculate the error using standard error propagation rules without 
taking care of the correlations between the ensembles [62j , we get errors that seem very low. 
I also repeated the FFS simulation only changing the length of the initial MD simulation. A 
decrease of 60% for the initial MD simulation resulted in a final result that is again 5 times 
smaller. For TIS, PPTIS, and RETIS the impact of this MD reduction would not even be 
noticed as it only effects the error in the flux term that is negligible compared to the error 
in the crossing probability. Fig. \7\ shows that the two different FFS crossing probabilities are 
similar at the start, but then start to deviate exponentially. The reason for this behavior is 
that the true set of reactive trajectories have an average kinetic energy distribution at the 
start that is considerably shifted compared to the equilibrium distribution. Therefore, a too 
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short MD simulation might not generate sufficient crossing points having a high velocity. As 
result, the FFS trajectories mainly climb up the barrier helped by the stochastic force instead 
of a high initial velocity. In Fig. [8] we compare five randomly selected crossing trajectories for 
TIS and for FFS. The FFS trajectories are clearly unrealistic as they are not symmetric in the 
(r, v) plane, which should be the case for a symmetric barrier. The velocities at the start are 
much lower than at the other side of the barrier. Moreover, two of the five trajectories, that 
were randomly picked from the 166 successful ones, start exactly from the same MD crossing 
point (from a total of 5260 crossing points). TIS, that is able to relax the history of the path, 
does not show this artefact. 

1.6, , , , , , , , , , , , 




-1.0 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1.0 1.2 

r 

Figure 8: TIS and FFS trajectories in the (r, v) plane 

These results have important consequences. It shows that attempts to use FFS for de- 
terministic MD [591l60j by applying some small level of stochastic noise can only work when 
inertia effects are not important. In other words, when the deterministic dynamics behaves 
effectively Brownian. This sampling problem is not unique to FFS, but to any splitting-type 
method such as weighted ensemble Brownian dynamics [721 [73], Russian Roulette [571155] . 
vector walking [71| , and S-PRES [75] , in which the equations of motion are only followed for- 
ward in time. There is not an easy solution for this problem. Still, this is very much desired 
for non-equilibrium systems for which there are no good alternatives. Possibly, the original 
nonequilibrium TPS approach of Crooks and Chandler ^76j could do better in this situation as 
it continues rebuilding trajectories from the beginning. Still, the set of starting points follows 
from a straightforward MD trajectory. Therefore, also this approach is likely to miss these rare 
initial points that have high potential to become reactive. Timereversal moves might alleviate 
the problem [lOj, but can only be applied for time-reversible dynamics and velocity-symmetric 
steady state distributions. The other possibility is to adapt the RC, X{x), to the phasepoint 
committor, for instance via an on-the-fly optimization scheme [66]. Such an approach would 
need to be carried out in full phasespace, which is presently not the standard. 
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However, redefining tlie intermediate interfaces to lie on the isocommittor surfaces alone 
will not be sufficient. If we keep the r-dependent definition for stable state A, the committor 
probability can jump from zero to one is a single timestep (at the point when it leaves state A 
with a high velocity). Hence, also the state definition A should be redefined in phasespace. For 
this particular system it seems intuitive to use constant energy curves X{x) = l/2mv'^ + V{r) as 
RC. FFS will probably work in that case.. However, it is yet unclear if it is practically feasible to 
design appropriate RCs using full phasespace in more complex systems. Present algorithms |14j 
have always assumed that is sufficient to use configuration dependent committor functions. 

On the other hand, the TIS methods seems to work properly using configuration dependent 
RCs. Only to ensure the stability of state A it is sometimes convenient to let Aq be velocity 
dependent [8j (For this system, the friction coefficient is sufficiently high to neglect kineticly 
correlated recrossing) . On the contrary, it seems that TIS and its variations do not necessarily 
improve when the RC equals the true isocommittor, a hypothesis that was postulated in 
Ref. [Sj. If interfaces are places at constant energy curves, the trajectories will become much 
longer than in the present case. 

9 Conclusions 

I have reviewed some dynamical rare event simulation techniques. The RF method is likely 
the most efficient approach when studying low dimensional systems for which an appropriate 
RC can easily be found. The most efficient implementation of the RF approach to calculate 
the dynamical factor is probably the EPF algorithm that is considerably more efficient than 
the more common transmission coefficient calculation schemes. However, even with this more 
efficient EPF approach, the RF efficiency will decrease exponentially with barrier height and 
inverse temperature, if a proper RC can not be found [39]. The TPS reactive trajectory sam- 
pling does not require a RC. However, a definition of a RC is still needed in the TPS rate 
calculation algorithm, which has been improved by the TIS and RETIS methodologies. The 
TPS/TIS/RETIS efficiency only scales quadratically with barrier height and inverse temper- 
ature when using an "improper" RC [39]. The RC insensitivity of these methods gives them 
a strong advantage compared to RF methods in complex systems. Of these methods, RETIS 
is significantly faster than the other two. However, its implementation is somewhat more 
difficult than TIS. PPTIS (and the similar milestoning) is not an exact method as it assumes 
memory loss beyond a traveling distance between two interfaces. Using this approximation, 
PPTIS is able to reduce the required path length considerably, which can be important for 
diffusive barrier crossings. FFS does not require information on the phasepoint density and is, 
therefore, ideally suited to study nonequilibrium events. The advantageous implementation 
and its apparent efficiency have made FFS a very popular method for equilibrium systems 
as well. However, the numerical study, presented here, shows that FFS has certain pitfalls 
that have not yet been reported. It shows that the RC sensitivity of FFS is even more trou- 
blesome than it is for RF methods. Present simulation studies have almost always assumed 
that RC are functions of configuration space alone. My example shows that an appropriate 
RC for FFS needs to be defined phasespace while configurational space would be sufficient 
for the RF method and the equilibrium path sampling algorithms TPS/(RE)TIS. Still, there 
are presently no alternative methods that can treat nonequilibrium processes and do not have 
the same problem. The fact that FFS and other forward MC methods get so easily trapped 
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towards unfavorable reaction paths, by missing an important orthogonal coordinate or ve- 
locity in an early stage, requires the uppermost caution when applying these methods and 
interpreting their results. In addition, this article poses challenges for developing improved 
nonequilibrium path sampling methods that are either less sensitive to a chosen RC or are 
able to find an appropriate phasespace dependent RC on the fly. 
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10 List of abbreviations 



EPF 


Effective Positive Flux 


FES 


Forward Flux Sampling 


MC 


Monte Carlo 


MD 


Molecular Dynamics 


PPTIS 


Partial Path Transition Interface Sampling 


RC 


Reaction Coordinate 


RE 


Replica Exchange 


RETIS 


Replica Exchange Transition Interface Sampling 


RF 


Reactive Flux method 


TI 


Thermodynamic Integration 


TIS 


Transition Interface Sampling 


TPS 


Transition Path Sampling 


TS 


Transition State 


TST 


Transition State Theory 


US 


Umbrella sampling 
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